A Cross-Validation Approach to Approximate Basis 
Function Selection of the Stall Flutter Response of a 
Rectangular Wing in a Wind Tunnel 


Sunil L. Kukreja* 

NA SA Dryden Flight Research Center, Edwards, California 93523 

Gareth A. Vick 

The University of Sydney, NSW 2006, Australia 

Thomas Andrianne* , Norizham Abdul Razak * and Grigorios Dimitriadis * 

University of Liege , Liege fOOO , Belgium 


The stall flutter response of a rectangular wing in a low speed wind tunnel is modelled 
using a nonlinear difference equation description. Static and dynamic tests are used to 
select a suitable model structure and basis function. Bifurcation criteria such as the Hopf 
condition and vibration amplitude variation with airspeed were used to ensure the model 
was representative of experimentally measured stall flutter phenomena. Dynamic test data 
were used to estimate model parameters and estimate an approximate basis function. 

I. Introduction 

The existence of nonlinearities in aeroelastic systems is known to significantly affect their dynamic be- 
havior, often inducing oscillations that cannot be predicted by linear theory . * 1 * One such type of oscillation 
is known as stall flutter. Traditionally, aeroelastic phenomena are considered to be dynamic responses 
resulting from interaction of structural, inertial and aerodynamic forces. Such interactions can lead to the 
well-known classical linear flutter, a catastrophic self-excited oscillation. However, several other incidents 
can arise from low speed aeroelastic interactions, such as stall flutter and galloping , 4 which are due to the 
presence of an underlying nonlinearity in the aerodynamic forces. These events frequently lead to Limit 
Cycle Oscillations (LCO). LCOs also result from other sources, such as structural or material nonlinearities. 

Stall flutter is a phenomenon that occurs when flow separates from and reattaches to the surface of a 
wing in a cyclic manner. The separation can be partial or complete over the entire wing’s surface. This 
alternation between stalled and attached flow is known as dynamic stall, a phenomenon that has been the 
subject of numerous experimental and theoretical investigations (see for example McCroskey et al , 5 Ericsson 
and Reding 6 and Spentzos et al 7 ). Stall flutter is the result of coupling of the vibration characteristics of a 
flexible structure with dynamic stall. 

Historically, several simple models of stall flutter have been developed, such as those by Ericsson. 
However, practically such models are quasi-steady, i.e. they are derived from steady lift and pitching moment 
curves of a wing. In the present work, a straightforward mathematical model of stall flutter is developed from 
measurement of dynamic responses. Abdul Razak et al ’ 10 conducted a set of experiments on a rectangular 
wing in a low speed wind tunnel that underwent stall flutter in pitch. These experiments are used to identify 
a nonlinear discrete-time model of stall flutter response of a wing and estimate the fundamental parameters 
of its dynamics and static function. 

The model was estimated using a cross-validation approach. 5 Model development via cross-validation 
is concerned with assessing whether a given nominal model can reproduce data from a plant, collected after 
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some initial experiments to obtain estimation data. The model validation problem is really one of model 
invalidation since a given model can only be said to be not invalidated with the current evidence. Future 
evidence may invalidate the model. 

The motivation for this study is to investigate whether an approximate basis function can be parame- 
terized and estimated using the £ 2 estimator, thus, avoiding a nonlinear least-squares estimation and search 
paradigm. While this work provides an ad-hoc approach to basis function parameterization it is not a re- 
placement for the traditional nonlinear least-squares approach. This work offers a coarse estimation, which 
may be adequate in most applications or provide a possible seeding for nonlinear search techniques such as 
the well-known Levenberg-Marquardt algorithm. ’ ’ 

II. Experiment Description 

Below we provide details of the experimental testbed and procedures used for data collection and analysis 
to obtain a basis function class. 

A. Support System 

The experiments ’ were conducted in the Multi-Disciplinary Low Speed Wind Tunnel of the University of 
Liege. The wind tunnel features two working sections, one for aeronautical and automotive applications and 
one for wind engineering applications. Experiments were carried out in the aeronautical working section, 
which has dimensions of 2m x 1.5m x 5m (width x height x length) and is capable of achieving airspeeds 
of up to 60 m/s in the tunnel’s closed loop configuration. The turbulence level was 0.15% on average at 
lower speed ranges. For stall flutter experiments a special support framework was assembled, which was 
used to fix extension springs. Ends of the springs were fixed to arms that acted as adaptors to the wing 
spar. The spar was connected to the center of the arm. Free-play was avoided by securely attaching the spar 
to the adaptor arms. The suspension springs allowed movement in all directions. The support mechanism 
described was designed to allow frictionless wing vibrations with known structural stiffness values. Eight to 
sixteen springs can be installed to suspend a wing and keep it in tension. Even numbers of springs were used 
to ensure symmetry in stiffness. The complete support system is shown in Fig. 1. 



(a) Drawing (b) Photograph in wind tunnel 


Figure 1. Wing and support design and construction 


B. Wing 

A rigid straight rectangular wing was constructed from aluminum and balsa wood. The wing chord was 
36.0 cm long and span 100 cm, resulting in an aspect ratio of 2.78. The skin consisted of a 0.5 mm thick 
aluminum sheet that was wrapped around balsa ribs. There were six ribs inside the wing, positioned 18.0 
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cm apart. The 2.0 cm thick ribs were held together by an aluminum rod acting as a main spar, located at 
37% of the chord. The airfoil selected for this study was the NACA 0018. This thickness ratio was chosen 
for two reasons. Firstly, it was high enough to allow the installation of instrumentation inside the wing. 
Secondly, it was thicker than the NACA 0012 airfoil investigated by Dimitriadis and Li. It was hoped 
that the increased thickness would give rise to different stall phenomena. This design allowed ballast to be 
attached internally. Mass of the wing including main spar, adaptor arms and instrumentation was 6.47 kg 
and its spring stiffness was 755 N/m. 

To measure the unsteady pressure distribution, 16 pressure tappings were drilled on the wing surface, at 
the mid-span point. These tappings were connected to 16 piezoresist ive pressure transducers 1 by means of 
plastic tubes. 1 They were positioned symmetrically on the wing’s upper and lower surfaces, as shown in 
Fig. 2. The phase lags associated with the tubes were discussed in. 



Figure 2. Position of pressure sensors, spar and center of gravity relative to the wing section 


Wing motion was measured using four accelerometers attached on the spring adaptor arms. Their 
frequency measurement range was between 0-300 Hz with accuracy of ±5% and their measurement range is 
±50 g. The pressure sensors and accelerometers were sampled at 1.0 kHz using a National Instruments PXI 
1010 data acquisition module managed by Labview 9 software. From these measurements, pitch and plunge 
displacements of the wing were reconstructed by numerically integrating acceleration values twice. There 
were no end- plates attached to the wing tips to allow good quality flow visualization could be obtained. As 
a result, flow over the wing was three-dimensional. 


C. Procedure 

Two parameters varied during these experiments are as follows: 

• Wind tunnel airspeed, V (from 0 m/s to a maximum of 26 m/s) 

• Wind-off static pitch angle, ao (from 11° to 16°) 

Experimental procedure consisted of setting the wind-off static pitch angle then gradually increasing 
airspeed. At each new airspeed, when wing response was stabilised, 4 seconds of acceleration and pressure 
response data was recorded. No external excitation was applied besides the natural turbulence of the wind 
tunnel. 

Before starting the wind-on tests, wind-off instrumented hammer tests were performed in order to de- 
termine the natural frequencies of the degrees of freedom of the wing. The frequency values for the plunge, 
pitch and roll degrees of freedom are tabulated in Table 1. The yaw, lateral translation and longitudinal 
degrees of freedom have higher natural frequencies. 

D. Data Analysis 

Accelerometer data were first analysed to obtain pitch, plunge and roll accelerations. These responses were 
band-pass filtered between 2Hz and 80Hz, using a rectangular frequency domain filter. They were then 
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Mode 

Frequency 

Damping ratio 

Plunge 

4.57Hz 

0.0033 

Pitch 

5.95Hz 

0.0087 

Roll 

13.10Hz 

0.0418 


Table 1. Natural frequencies and dampings at U = 0 m/s 


integrated in time twice to estimate the pitch, plunge and roll velocities and displacements. Integration was 
carried out in the time domain using a central difference scheme. 

Pressure data were integrated around the wing’s central section to obtain sectional aerodynamic lift 
and moment coefficients. These signals were band-pass filtered between 2Hz and 80Hz. Note that lift and 
moment coefficients are sectional and, therefore, do not represent the tota] lift and moment acting on the 
wing. In this work, the sectional moment coefficient is used to choose a basis function for identification of 
the wing’s response. This choice is based on variation of sectional lift and moment coefficient values with 
pitch displacement and velocity responses. 


Moment coefficient around flexural axis 



Figure 3. Static pitching moment variation with angle of attack at two different Reynolds numbers. 


III. Results 

Only a subset of the experimental results is presented in this work. For a more complete description see 
Abdul Razak et al. 0 Static aerodynamic pitching moment results are presented first, followed by dynamic 
responses. 

A. Static Results 

Figure 3 shows the total static aerodynamic pitching moment coefficient variation around the pitching axis 
with pitch angle. The data shown in the figure were obtained using a force and moment balance setup. The 
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wing was placed at different pitch angles and its total lift and drag were measured using six force gauges; 
aerodynamic moment around the pitching axis was measured using a single torque gauge. Static experiments 
were repeated at two Reynolds numbers representative of subsequent dynamic experiments. The pitching 
moment coefficient is defined as 

711 

m 1/2 pV 2 Sc 1 j 

where m is the pitching moment, p is the air density, V is the wind tunnel airspeed, S the wing surface area 
and c is the chord of the wing. 

Figure 3 illustrates that the moment curve is linear up to a pitch angle of approximately 14°. Nose up 
pitching moments are defined as positive and, hence, the effect of the pitching moment is destabilising. At 
pitch angles higher than the linear limit the pitching moment becomes less destabilising. 

The static pitching moment coefficient data can be curve-fitted using a function of the form 

c m = hiaexp(-h 2 a h3 ) (2) 

where c m is the non-dimensional aerodynamic pitching moment coefficient and a is the pitch angle. Note 
that Fig. 3 plots a curve-fit based on Eqn. 2, where hi = 0.4011, = 707.5641, h% = 6 and a is in rad. The 

curve- fit provides a good approximation up to pitch angles of 24° . 

Since pitch responses were obtained from integrating acceleration signals, the mean values of these re- 
sponses were lost. However, mean data can be reconstructed using the results in Fig. 3. It is assumed that 
the mean value of pitch response is equal to the static deformation of the wing’s pitch angle at the same 
airspeed. This deformation is given by the static equilibrium equation 

kcx{oL e — C^o) = 2^ ^ C Crn static( a e ) (3) 

where k a is wing stiffness in the pitch degree of freedom, a e is mean pitch angle at airspeed V and Oig is 
wind-off pitch angle. This equation is nonlinear but can be solved iteratively or graphically. 

B. Dynamic Results 

Dynamic tests were performed after the static tests, replacing the force balance setup with the spring 
assembly shown in Fig. 1. For all selected wind-off static pitch angles, at low airspeeds the response was 
noisy and of low amplitude. For a static pitch angle of 11°, significant high amplitude vibrations occurred 
at airspeeds higher than 25 m/s. This is a classical flutter phenomenon. The highest airspeed tested was 
limited to 25.5 m/s for safety reasons. 

For the other static pitch angles, lower amplitude vibrations were encountered at lower airspeeds. Vi- 
brations occurred primarily in pitch degree of freedom, with very little motion in other degrees of freedom. 
The amplitude of these vibrations increased gradually with airspeed, particularly for the highest static pitch 
angles of 14° — 16°. Figure 4 shows variation of pitch angle response amplitude with airspeed for all tested 
angles of attack. 


IV. Discrete-Time Nonlinear Model Representation 

In this work, the construction of identified models describing the stall flutter response of a wing is 
considered. Therefore, low amplitude responses to turbulence at low airspeeds and the classical flutter 
phenomenon at ckq = h° is ignored. As mentioned previously, stall flutter responses predominantly involve 
pitch degree of freedom. Therefore, the phenomenon can be modelled as a single degree of freedom vibration. 
The equation of motion for the wing can be written as 

I a a(t) + c a d(t) + k a a(t) = ^ pV 2 Sc c TO (a(U), a(ti), d(ii)) (4) 

where t is time, a is pitch angle, I a is wing’s moment of inertia around the pitch axis, c a is wing structural 
damping in the pitch degree of freedom, p is air density and c is wing chord length. In Eqn. 4, the pitching 
moment coefficient is expressed as a nonlinear function of the current and past values of system states. The 
motion started at time t = 0, and t\ takes values from 0 to t. Dependence on past values of system states 
reflects the well-known property of infinite dimensionality of unsteady aerodynamics. 
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Figure 4. Pitch response amplitudes at all angles of attack and airspeeds 


Equation 4 can be simplified if it is assumed that the added mass effect is small compared to the moment 
of inertia of the wing’s structure, even when flow is stalled. In this case, Eqn. 4 becomes 

l a a(t.) + c a a(t) + k a a(t) = ^ pV 2 Sc c m (a(ti), a(ii)) (5) 

To reflect experimental conditions in the wind tunnel, an additional excitation force term can be added 
to model the effect of free stream turbulence. If it is assumed that this effect is much smaller than the 
self-excited stall flutter oscillations, the equation of motion can be written as 

I a a(t) + c a a(t) + k a a(t,) = ^ pV 2 Sc c m (a{ti), a{ti)) + f(t) (6) 

where f(t) is a random excitation force of small amplitude. 

A. Nonlinear Structure of Moment Coefficient 

Using the present approach, the moment coefficient is an unknown nonlinear function of pitch displacement 
and velocity at all times since the onset of motion. To choose an approximate basis and model form for this 
nonlinear system, sectional pitching moment coefficient signals measured using the pressure tappings were 
utilized. These signals do not represent a complete 3D aerodynamic pitching moment acting on the wing. 
However, it is assumed that the sectional moment coefficient has the same nonlinear form as the 3D moment 
coefficient and, therefore, can be used to provide a good approximation. 

Figure 5(a) shows the time response of a, a and c m for ao = 14° and V = 11.5 m/s. The oscillation 
amplitudes are very small (less than 1° in pitch displacement) and the signals feature significant randomness. 
Figure 5(a) depicts a typical response to the wind tunnel’s free stream turbulence, in the absence of any self- 
excited phenomena. Figure 5(b) displays c m on the 2 -axis against a and a on the x and y axes, respectively. 
This plot is very messy and without any clear structure. 

By contrast, Fig. 6 shows the time response and phase plane plots at the same wind-off pitch angle 
but a higher airspeed, V = 20.1 m/s. In this case, the signals are significantly more ordered and have 
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(b) Phase plane 

Figure 5. Time responses and phase plane plots of a , d and c m at low airspeed 


a single principal frequency component. Nevertheless, the moment coefficient varies more stochastically 
than pitch displacement and velocity responses. This phenomenon is due to the fact that the separation 
and reattachment of the boundary layer are quasi-repeat able processes. During each cycle, separation and 
reattachment occur at slightly different phase values. As a consequence, the phase plot of Fig. 6(b) is noisy 
but has a clear structure. Furthermore, it is repeatable within some confidence bounds. 


V=20.1m/s 




(a) Time responses (b) Phase plane 

Figure 6. Time responses and phase plane plots of a, a and c m at high airspeed 

The objective of this section is to develop a discrete-time model representation for c m , which is capable of 
representing phase-plane plots such as Fig. 6(b). If no nonlinearity exists or the system is perturbed around 
a sufficiently small operating point, an ARMA (Auto- Regressive Moving Average) model structure would 
suffice and is represented as, 


c m (n) = 70 a(n) + 71 a(n - 1 ) + 72 a(n - 2 ) + . . . + e(n) + S 1 e(n — 1 ) + S 2 e(n - 2) . . . ( 7 ) 

where n is a discrete-time instance, e(n) is a noise sequence and 7 y, Si are coefficients to be estimated. 
However, such a model cannot represent the correct geometric shape of the manifold in Fig. 6(b) and will 
result in a linear model for the entire wing that is incapable of predicting LCOs. Therefore, a nonlinear 
model must be formulated. 

To choose an appropriate nonlinear model structure for the manifold of Fig. 6(b), it is observed that the 
plot is a 3 D curve. However, if the manifold is rotated appropriately, it can be seen that it lies on a curved 
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V=20.1m/s 




Figure 7. Figure 6(b) plotted at different azimuth and elevation angles 


2D plane . Figure 7 shows the same manifold at two different sets of azimuth and elevation angles. This 
illustrates that for two combinations of azimuth and elevation, the manifold is an approximately 2D curve. 
This phenomenon is common to all measured LCOs, although rotation angles differ for each case. Notice 
that the plot in Fig. 7(a) resembles the static pitching moment curve in Fig. 3 at pitch angles between 10° 
and 20°. Therefore, to obtain an estimate of c m , the cycle a(a) can be mapped onto a curved plane, as 
shown in Fig. 8. The shape of the plane is inspired by the function used to curve fit the static pitching 
moment coefficient data, as shown in Eqn. 2. The coordinates of that plane are given by 


x = RR Z 


a/w x i 

dijWy 

— a exp (— a A /w w )/w z 


+ bw 


(8) 


where x = [x y z] T are coordinates of the mapping plane in a-a-c m space, w x , w y and w z are normalisation 
coefficients, 6 is a displacement coefficient and v is a unit vector perpendicular to the direction of the curvature 
of the plane, given by 

v = [cos <j) sirup Of (9) 

where <p is the azimuth angle that collapses the c m manifold into a 2D curve, in this case —40° as in Fig. 7(a). 
K z is the azimuthal rotation matrix given by 


r. 


/ cos <j> — sin^> 0 

sin (j) cos <f> 0 

\ 0 0 1 


(10) 


and R rotates the mapping plane around an axis aligned with v by an angle equal to the elevation angle ^ 
required to collapse the c m manifold to two dimensions, in this case —4°. R is given by 


R = 


^ cos^ + cos 2 <p 2 {l — cos^) 
cos (f) sin 4>{ 1 — cos 'ip) 
y — sin <p sin pj 


cos <j> sin <f>(l — cos^p) 
cos %p + sin 2 <f>( 1 — cos r ip) 
cos <f> sin y? 


sin 4> sin ^ ^ 
cos (f) sin ip 
cos^ / 


(ii) 


From Eqn. 8, it is clear that a good approximation of the geometric plane on which e m lies has been 
obtained by means of a simple nonlinear function. Note that in the case of Eqn. 8 the exponent of a inside 
the exponential is 4 as opposed to 6 in the static pitching moment case. In the sequel, a generalized form of 
this basis function is used to develop a robust model description. 
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0 



Figure 8. Moment coefficient manifold and mapping plane 


B. Constraints 


The nonlinear function c m must satisfy some known constraints. First, the system undergoes a Hopf Bi- 
furcation at a critical airspeed, V c , such that a periodic solution branch appears in the system’s bifurcation 
diagram at that airspeed. For a single degree of freedom system of the form of Eqn. 5 or 6, the Hopf bi- 
furcation can be described by considering stability of the equilibrium solution a = a e . The stability of this 
solution is identical to the stability of the linear system 


J a d(t) + 


c a a(t) + k a a(t) = -pV 2 Sc 



a T 


a e ,V 


dCm 

da 



(12) 


For a Hopf bifurcation to occur, total damping of the linear system above must be equal to zero. This yields 
a bifurcation condition for the form 


OCm 

da 


2 c 0 


a e ,V c 


pV c 2 Sc 


(13) 


Second, the nonlinearity must satisfy the limit cycle amplitude and fundamental frequency. At airspeeds 
higher than the bifurcation airspeed the system undergoes LCOs whose amplitude and frequency have been 
measured. Therefore, the nonlinearity must be capable of matching these measured values. To describe this 
condition, a first order Harmonic Balance expansion is applied to the equation of motion. It is assumed that 
the LCO response of the nonlinear equation of motion can be written as 


a (t) = A sin cut (14) 

where A is LCO amplitude and cc is frequency. Substituting this solution into the equation of motion yields 


— Aou 2 I a sin out + Aouc a cos out + Ak a sin out = -pV 2 Sc (ao + cq cos out + b\ sin out) 
where ao, a% and b\ are coefficients of a Fourier series given by 

f*27T / C xJ 


( 15 ) 


a 0 


f- Z7T j UJ 

I Cm 

Jo 


A sin out, A cos out ) dt 
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^ pzir/co 

Cbl = — Cm (Asiuout, A cos out) COS outfit 

K Jo 

j'2'K I u) 

bi = — / c m (A sin out ^ A cos out) sin outdt. 

n Jo 


Substituting into Eqn. 15 and carrying out harmonic balancing yields 


/»27T / U> 

/ Cm (A sin out, A cos out) dt = 

Jo 

f‘ C l'K j U) 

/ Cm {A sin out, A cos out) sin outdt = 

Jo 

/•27T fu> 

/ c m (. A sin Ccff, A cos out ) cos outdt = 

Jo 


2io(x e 

pV 2 Scou 

2ttA 

p\ /2 Scou 

2ttA 

pV 2 Scou 


{—^ 2 da + fca) 


(16) 


Conditions 13 and 16 can be used in order to select a suitable set of basis functions for the nonlinear term 

Cm- 


C. Choice of Nonlinear Function 

Satisfying the constraints above yields a suitable form of nonlinear function c m represented as 

OL 

c m =piaexp(-p 2 a P3 ) +P4-. (17) 

The linear term in a is necessary to satisfy condition 13. Nevertheless, this model is restrictive because the 
exponent p 3 must be an even integer, if the function is to remain real and finite at negative values of a. A 
generalization of the model, that allows odd and non-integer exponents can be written as 

OL 

Cm =P\OL exp (— p 2 K 3 |) +P±y (18) 


where | | denotes the absolute value. 


D. Full Model 


Using the formulation of Eqn. 18, the equation of motion of the wing in pitch becomes 


IaA{t) + c a a(t) + k a (ct(t) 


a 0 ) = ^pV 2 Sc p 1 aexp(-p 2 |a P3 |) + ^ pV 2 Sc 


(19) 


where the centering of the pitch spring around the wind- off pitch angle is taken into account. 

Equation 19 can be posed in discrete- time form using a central difference approximation, centered around 
time index n — 1. The pitch acceleration and velocity terms are written as 


a 

a 


a(n) — 2 a(n — 1) + a(n — 2) 
At 2 

a(n) — a(n — 2) 

2A t 


where the notation a(n) denotes a(nAt) and At is the time step between two successive experimental 
measurements, in this case equal to 0.001s. The discrete-time form of the equation of motion is then 


f k_x C a 

VAf2 


1/2 pVScp 4 \ , 
— CL[n 

2 At J v ' 


k a a 0 + 



a(n — 1) 


+ I “AC + 


c a — 1/2 pV Sc p 4 
2At 


a(n — 2) + -pV 2 Sc pi a(n — 1) exp (— p 2 \a(n — 1) P3 1) . 


( 20 ) 
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Wind- off dynamic tests provide values of several coefficients in this formulation, namely, f a . c a and k a . 

Furthermore, p 4 is known from Eqn. 13, i.e. 


Pa 


2c a 

pV c Sc 


(21) 


The exponent p% is not known a priori. Unfortunately, integrals in the first and second conditions of Eqn. 16 
cannot be evaluated analytically However, they can be evaluated numerically such that 


2ir /cj A t 

y pia(n) exp (— p 2 | a(n) P3 1) At 

n=0 

2n/ujAt 

y pia(n) exp (— p 2 \a(n) P3 1) sin (urn At) At 

n = 0 


2na e 

pV 2 Scuj 


2nA 

pV 2 Scuj 


( —CJ 2 I a + ka) . 


(22) 


Therefore, a complete model of the wing can be obtained by estimating the unknown parameters of Eqn. 20 
given the constraints 22. 


V. Best Fit Model via Cross-Validation 

Although Eqn. 20 provides a compact model description, it yields a non-trivial estimation problem. 
A typical estimation procedure would involve a nonlinear least-squares approach because of the unknown 
coefficients associated with the exponential term. In this paper, we chose to parameterize a set of basis 
functions using a priori knowledge from the present dynamic experiments and using a cross-validation 
technique to select a best fit model which contains an approximate basis function. 

A. Predicitive Capability Using an ARMA Model Formulation 

In many cases although a nonlinear model is suspected or assumed from finite element modelling efforts, 
experiments are conducted about a small operating region due to safety concerns and hardware limitations. 
This often renders the nonlinear effect unobservable and, hence, unidentifiable. 21 As such, we first estimate 
a linear model, namely Eqn. 20 with the basis function truncated to assess whether a linear model is capable 
of predicting the observed output before proceeding to postulate a more advanced nonlinear one. 

Specifically we used one of eleven datasets for estimation and the remaining ten to validate the linear 
model’s predictive capability. The accuracy of a Ustep-ahead predictor was quantified by computing the 
%Fit as 

%Fit = (1- X 100 0 3 ) 

V N 2n= 1 ( a ™) / 

where a is the predicted output. Figure 9 shows a typical cross-validation output for the postulated ARMA 
model. The figure illustrates that a linear model is not capable of accurately reproducing the measured data 
with high accuracy. For all cross-validation datasets the ARMA model was only able to predict approximately 
75% of the observed output variance. This result provides justification to assess whether a nonlinear model 
can provide improved predictive capability. 

B. Predicitive Capability Using a Nonlinear Model Formulation 

We evaluated the ability of the nonlinear model given by Eqn. 20 to accurately predict the observed output. 
The parameterization chosen for the exponential terms, P 2 and ps is given in Table 2. The values of p 2 and 
ps were selected by trial and error as ones that provided good curve-fits with a maximum error of less than 
10%. For estimation and model selection, we again used one output as estimation data and the remaining ten 
as validation data. A model (Eqn. 20) was estimated for each parameterization giving 16 models for every 
validation data record. Each model’s output was computed using a Avstep-ahead predictor and accuracy 
quantified as %Fit (Eqn. 23). The model that yielded the highest %Fit was deemed the best fit model. This 
procedure was repeated for each validation output. 
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Typical Cross-Validation Without Nonlinear Term 
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Figure 9. Cross-Validation: Typical linear model prediction (dash-dot line superimposed on top of the 

measured output (solid-line 


P2 

0.8 

0.9 

1.0 

i.i 

1.2 

1.3 

1.4 

1.5 

1.6 

1.7 

1.8 

1.9 

2.0 

2.1 

2.2 

2.3 

P3 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

1.1 

1.2 

1.3 

1.4 

1.5 

1.6 


Table 2. Values of p 2 and ps used to parameterize the basis function in Eqn. 20. 


The results produced the same best fit model for each validation data indicating the estimated model is 
robust. The selected model was consistently parameterized with the values, p 2 =2.0 and ps = 1.3. 

Figure 10 provides a typical cross-validation output for the nonlinear model. The figure illustrates that 
this nonlinear model is capable of accurately reproducing the measured data with high accuracy. For all 
cross-validation datasets, the discrete-time nonlinear model was able to predict 99.2% of the observed output 
variance. 

Figure 11 summarizes the cross-validation %Fit’s for each validation data, The results show that the 
predicted output, for these parameter estimates, account for a large portion of the variance. This indicates 
that the model parameters explain the measured data well. 

VI. Discussion 

Our results show that a linear model is not capable of reproducing the observed output for this model of 
stall flutter response. In many cases, researchers have not been able to quantify nonlinearities in their data, 
possibly due to it being collected about a small operating point, making the nonlinear effect unobservable. 

Although our model is nonlinear, we posed a linear statistical model to implement a best fit model 
selection approach via cross-validation data to avoid using a nonlinear least-squares estimation approach. 
This approach was feasible due to having a priori insight to good parameterization (see Tab. 2) to be used 
with this model. Although this methodology provided a high percent fit for all validation sets, since the 
underlying model is nonlinear and we used a coarse parameterization, the estimated parameters may be at 
a local minimum. Nevertheless, the estimated model could be used in applications where the concern is to 
achieve high fidelity simulations or as a seeding in a nonlinear least-squares framework possibly facilitating 
rapid convergence to a global minimum. 

In this study, we only considered one class of an infinite number of basis functions. It may be possible 
to reformulate this problem as a proper linear statistical model, thus, avoiding the need to have insight to 
“good” parameterization of the exponential term or other similar basis functions. In addition, recasting this 
problem in a linear framework offers the advantage of not converging to a local minimum. The results of 
this study will be presented in a future paper. 
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Figure 11. Cross-validation: %Fit versus experimental trial. 


Cross-Validation Fit 


2 4 6 8 10 

Experiment Number 


VII. Conclusion 

This study illustrates that nonlinear effects are present in the staff flutter response of an aeroelastic 
wing structure and that a nonlinear model is capable of modelling the dynamics with high accuracy. We 
showed that this nonlinear model can be accurately posed as a linear statistical model allowing i 2 estimation 
techniques to be used and, thus, avoiding a full nonlinear least-squares solution. These results may have 
utility for initialization of nonlinear least- square problems providing a global minimum in fewer steps than 
traditional approaches. 
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